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1.  Introduction 


Rapid  characterization  of  complex  urban  environments  via  physics-based  numerical  modeling 
will  likely  provide  important  infonnation  to  U.S.  Army  Soldiers  on  the  performance  of  advanced 
sensors,  as  well  as  the  effectiveness  of  computer  aids  to  increase  situational  awareness. 

However,  two  current  (and  extensive)  surveys  of  the  literature  (1,2)  indicate  that  computer 
simulation  of  wind  flow  and  temperatures  around  complex  urban  structures  have  most  often  been 
achieved  via  intricate  computational  fluid  dynamics  (CFD)  codes,  which  are  (as  a  rule)  quite 
computationally  intensive.  For  example,  CFD  codes  can  require  1  to  8  hours  of  execution  time 
on  multiprocessor  super-computers.  In  addition,  many  of  the  CFD  models  described  in  current 
literature  have  been  in  use  and/or  in  development  for  10  or  more  years.  Also,  CFD  models  are 
generally  cumbersome  to  modify  and  debug.  Hence,  something  inbetween  is  needed,  e.g., 
something  that  is  more  computationally  efficient  but  has  enough  flexibility  to  apply  to  the  types 
of  field  tests  that  are  envisioned  for  future  efforts.  Nevertheless,  it  may  be  useful  to  investigate 
CFD  model  frameworks  to  gain  a  better  understanding  of  the  considerable  task  that  is  undertaken 
when  attempting  to  develop  such  software.  Then,  one  can  begin  to  explore  alternate  model 
frameworks,  which  are  reliable,  rapid,  and  robust  to  simulate  meteorology  and  turbulence  in 
urban  environments  (e.g.,  to  study  atmospheric  effects  on  acoustic  propagation  or  free-space 
optical  communications). 

This  paper  initiates  the  investigative  process  by  presenting  a  critical  assessment  of  six  CFD 
model  frameworks.  Section  2  gives  a  brief  review  of  the  different  model  simulations  of  wind 
flow  and  the  thermal  microclimate  around  single  and/or  multiple  buildings.  Two  of  the  selected 
models  account  for  one  or  more  embedded  tree  arrays.  Additional  information  (if  available)  is 
provided  regarding  numerical  methods,  physics/turbulence  algorithms,  model  domain,  grid 
spacing,  time-step,  runtime,  and  computer  platfonn.  Section  3  discusses  some  of  the  difficulties 
and/or  deficiencies  with  each  modeling  approach.  Section  4  gives  a  summary  and  conclusions. 
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2.  Model  Survey 


The  numerical  models  described  in  this  section  were  selected  via  an  electronic  internet  search  of 
the  most  current  literature.  Selections  were  based  on  accessibility  of  information  regarding 
model  type,  numerical  methods,  physics/turbulence  algorithms,  grid  spacing,  time  step,  model 
domain,  runtime,  and  computer  platfonn.  Table  1  presents  a  comparison  chart  for  the  six  CFD 
models,  which  contains  these  kinds  of  data.  The  data  in  table  1  show  that  the  selected  CFD 
models  make  use  of  different  numerical  methods  (e.g.,  finite  difference,  finite  volume,  and  finite 
element)  to  mathematically  solve  the  equation  set.  Also,  the  selected  CFD  models  employ 
combinations  of  different  physics/turbulence  approaches  [e.g.,  Reynolds  Averaged  Navier- 
Stokes  (RANS),  large  eddy  simulation  (LES),  and  kinetic  energy — dissipation  ( k  -  s  )  turbulence 
models]  to  resolve  the  computed  fields. 
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Tablel.  Comparison  of  CFD  models. 


Computational  method/ 
Physics  approach 

Author(s) 

Features 

Turbulence 

Grid  spacing  /  Time  step 

Runtime 

3D  Finite  Difference 
(Incompressible  flow) 

Bruse  and  Fleer  (1998) 

Multiple  buildings; 
Embedded  array  of  trees 

k  -  e  model 

1.5  order  closure 

61  x  56  x  25  grids 

Ax  =  Ay  =  5m;  A z  =  4m;  At  =  10s 

6.0  hrs 

3D  Finite  Difference  -  RANS 
(Pseudo-compressible  flow) 

Wang  et  al.,  (2004) 

Pedestrian  winds  around  tall 
buildings 

k  -  s  model 

5 1  x  1 63  x  7 1  grids 

200  m  x  648  m  x  280  m 

Ax  =  Ay  =  Az  =  4m;  At  =  0.04s 

(t  =  20  min) 

(~8.3  hrs) 

3D  Finite  Control  Volume 
(Incompressible  flow) 

Paterson  and  Apelt  (1989) 

Single  prismatic  building 

k  -  £  model 

Non-uniform  staggered  grid 
Steady  state 

15  min 

(IBM  3083E) 

Johnson  et  al.,  (1997) 
Flerbert  et  al.,  (1998) 

Urban  canyon  winds  and  thermal 
microclimate 
(two  buildings) 

k  -  £  model 

240  m  x  632  m  x  100  m 

1  -2  m  grid  inside  canyon 

15-20  m  grid  outside  canyon 

45  min. 

(8  processor  super 
computer) 

3D  Finite  Volume  -  RANS 
(Incompressible  flow) 

Kim  and  Baik  (2004) 

Baik  et  al.,  (2003) 

Multiple  building  array 

RNG  k-  e  model 

Non-uniform  staggered  grid 

101  x  101  x  41  cells 

63.1  m  x  63.1  m  x  28.5  m 

At  =  0.055 

t  =  20  min 
(~  6.7  hrs) 

3D  Finite  Volume  -  LES 
(Compressible  flow) 

Pullen  et  al.,  (2004); 

Boris  (2002) 

Multiple  buildings 
(Central  business  district) 

MILES  model 

860  x  580  x  40  grids 
(Washington  DC) 

360  x  360  x  55  grids 
(Chicago) 

Ax  =  Ay  =  Az  =  6m;  At  =  0.36s 

8.0  hrs 

(16  processor  super 
computer) 

3D  Finite  Element  -  RANS 
(Incompressible  flow) 

Calhoun  et  al.,  (2004) 

Single  complex  building; 
Nearby  array  of  trees 

similarity-^ 

model 

1.0-2. 5  x  106  grid  pts. 

400  m  x  400  m  x  80  m 
Non-uniform  grid 
(finest  grid  =  1  m) 

~  1.0  hr 

(128  processor  super 
computer) 

2.1  Three-Dimensional  (3-D)  Finite  Difference  Model  (Multiple  buildings  with  an 
embedded  array  of  trees) 

The  paper  given  by  Bruse  and  Fleer  (3)  describes  a  non-hydrostatic,  3-D,  microscale,  numerical 
model  (called  ENVI-met)  for  surface-plant-air  interactions  in  and  around  urban  structures.  The 
model  ENVI-met  solves  the  basic,  incompressible,  Navier-Stokes  equations  forward  in  time  via 
finite  difference  numerical  methods.  The  main  model  equations  as  presented  by  Bruse  and  Fleer 
(3)  are  as  follows: 

For  i  =  1,2,3 


du  du 

- 1-  U ;  — 

dt  dXj 


-®  +  Kn 
dx 


dxf 


+ 


f(v-v  )-Su 


(1) 


dv  dv 

- 1"  U  :  - 

dt  dXj 


-SJ-  +  K, 

dy 


f  rs2  ^ 
O  V 


dxf 
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dw  dw 

- h  U  :  - 

dt  dxt 


dp' 

dz 


+  K„ 


f  ' 
o  w 

dxf 


+  g 


#(z) 

°refiZ) 


-3„ 


(3) 


du  dv  dw 

- 1 - 1 - 

dx  dy  dz 


(4) 


dG 


dG 


- 1"  U  :  -  —  K 


f  d2G^ 


dt 


dx: 


dxf 


+  Qh 


(5) 


and 


da  da 

—  +  U:  —L 

dt  dx, 


=  K, 


d2q 

dxf 


+  Qq 


l  J 


(6) 


Here,  t  is  the  independent  variable  time,  u,  v,  and  w  are  the  mean  wind  velocity  components, 
is  the  i-component  of  the  wind  velocity  vector  (ii\  =  u,  a 2  =  v,  M3  =  w),  and  x,  is  the  i-component 
of  the  position  vector  (jci  =x,X2=  y,  X3  =  z).  In  addition,  p  ’  is  the  local  pressure  perturbation,  Km, 
Kh,  and  Kq  are  the  turbulent  exchange  coefficients  for  momentum,  heat,  and  moisture, 
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respectively,/ (=  104  sec1)  is  the  Coriolis  parameter,  ug  and  vg  are  the  geostrophic  wind 
components,  6  is  potential  temperature,  and  q  is  specific  humidity.  In  equations  1  through  3,  the 
source/sink  terms  (Su,  Sv,  and  Sw)  describe  the  loss  of  wind  speed  due  to  drag  forces  from 
vegetation.  In  equations  5  and  6,  Qh  and  Qq  are  the  source/sink  terms  for  atmospheric  heat  and 
moisture,  respectively.  Two  additional  equations  describe  the  1 .5  order  closure,  k  -  s 
turbulence  sub-model  {4,5).  They  are  as  follows:  For  i  =  1,2,3 


8E 


dE 


- r  u  i  —  —  K  r 


dt 


8X: 


( d1E^ 

&,2 


+  Pr-  Th  +  Qe  -  s  , 


(V) 


and 


8s 


8s 


- h  U  ;  -  —  K 


k  o2  'N 
8  s 


dt 


8X: 


8x{ 


+  ci  c3^Th-c2^-  +  Qe 


(8) 


Here,  E  is  the  local  turbulence  (i.e.,  turbulent  kinetic  energy  (t.k.e.),  where  E  =  k  =  m,m,  /2  ),  s  is 
the  t.k.e.  dissipation  rate,  Pr  is  the  mechanical  production  of  t.k.e.,  Th  is  the  buoyancy 
production  of  t.k.e.,  KE  and  AT  arc  exchange  coefficients,  QE  and  Qs  are  source/sink  terms,  and  c;, 
C2,  and  cj,  are  numerical  constants. 

To  solve  the  combined  advection-diffusion  equations,  the  alternating  directions  implicit  method 
(ADI)  and  an  upstream  advection  scheme  is  used.  Dynamic  pressure  is  removed  from  the 
equations  of  motion  and  calculated  separately  from  the  Poisson  equation.  The  model  can 
simulate  wind  field  modifications  around  solid  boundaries  like  walls  as  well  as  modifications 
through  porous  media  like  trees.  The  ENVI-met  model  contains  sub-models  for  the  mean  wind 
flow,  temperature  and  humidity,  turbulence  and  kinetic  energy  processes,  radiative  fluxes,  soil 
and  vegetation  interactions,  and  ground  surface  and  wall(s)  interactions.  Bruse  and  Fleer  (3) 
provide  an  interesting  case  study  to  show  changes  in  local  wind  flow  and  calculated  temperature 
fields  through  a  typical  central  business  district  (figures  1  through  3). 
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OZU  Asphalt  Road  O  New  Trees 

Z  Pavement  □  New  Parkarea 

K1  Building  (Height  in  m) 


Figure  1.  Schematic  of  the  building  geometry  for  an  ENVI-met 
model  case  study.  Outer  buildings  are  24  m  in 
height  and  center  buildings  are  15  m.  Some  trees 
are  planted  along  the  upper  street  canyon  (from 
Bruse  and  Fleer  [5]). 


6 


Figure  2.  An  ENVI-met  model  calculation  of  the  horizontal  wind  field  at  z  =  2  m,  where  the  initial  wind 
direction  is  9=  90°  (from  Bruse  and  Fleer  [5]). 
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Figure  3.  An  ENVI-met  model  calculation  of  the  (normalized)  temperature  field  at 
z  =  2  m.  The  grey  area  indicates  a  central  park  with  trees.  Different 
surface  temperatures  result  due  to  shading  by  buildings  and  trees.  These 
effects  produce  a  pattern  of  warmer  (w)  or  colder  (c)  areas  than  the  average 
reference  temperature  at  the  same  level  (from  Bruse  and  Fleer  [5]). 


The  main  time  step  for  the  ENVI-met  model  is  At  =  10  s.  Smaller  time  steps  are  used  for  the  k—s 
turbulence  model.  Even  grid  spacing  (Ax  =  Ay  =  5  and  Az  =  4m)  is  used  in  each  direction, 
however,  the  lowest  grid  cell  above  ground  is  split  into  five  cells  with  size  A zg  =  0.2Az  to 
increase  accuracy  in  calculating  surface  processes,  e.g.,  the  surface  radiation  and  energy  budget. 
In  the  case  study  described  above,  the  ENVI-met  model  contained  61  x  56  x  25  grid  points  (300 
x  275  x  96  m).  The  model  calculation  takes  approximately  6  hours  to  resolve  the  computed 
fields. 

2.2  Three-Dimensional  (3-D)  Finite  Difference  Model  -  RANS  (Pedestrian  winds  around 
tall  buildings) 

Wang  et  al.,  ( 6)  describe  a  3-D,  microscale,  wind  flow  model  [called  PUMA  (Peking  University 
Model  of  Atmospheric  Environment)]  to  calculate  pedestrian  winds  around  tall  buildings.  The 
PUMA  model  is  based  on  the  Reynolds  Averaged  Navier-Stokes  (RANS)  equations,  where  the 
atmosphere  is  assumed  to  be  neutral,  i.e.,  without  thermal  effects.  (Note  that  implementing 
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energy  equations,  heat  flux  equations,  source/sink  tenns,  and  buoyancy  effects  frequently 
demand  additional  computer  time  and  resources,  which  the  developer  may  not  consider 
necessary  to  achieve  acceptable  model  results).  The  main  model  equations  as  presented  by 
Wang  et  ah,  (6)  are  as  follows:  For  i,j  =  1,2,3 


Giij 

dxj 


=  0, 


(9) 


du:  8u ,  1  dp 

dt  7  dxj  p  dxj 


(10) 


where  u,-  is  the  mean  velocity  component  in  the  7-th  direction,  p  is  the  air  density,  and  p  is  the 
fluctuating  pressure.  The  Reynolds  stress,  /?(/  =  uluJ  ,  represents  the  effects  from  turbulence. 

Familiar  summation  notation  is  used  in  equations  9  and  10  (and  elsewhere  in  this  section  of  the 
report).  Tunick  (7)  provides  several  useful  examples  to  demonstrate  how  the  complete  equation 
set  can  be  extracted  by  expanding  the  summation  indices. 


The  wind  flow  equation  is  integrated  forward  in  time  via  finite  differencing,  although  the  virtual 
compress  (pseudo-compressible  flow)  method  as  described  by  Chorin  (5)  is  adopted  to  solve  for 
the  pressure  field.  The  Reynolds  stress  (turbulence)  is  solved  via  a  modified  k  -  £  turbulence 
model  as  described  by  Jones  and  Lauder  (9).  The  equations  for  the  t.k.e.  ( k )  and  its  dissipation 
(s)  as  presented  by  Wang  et  ah,  (6)  are  as  follows:  Using  regular  summation  notation  for  i,j  = 
1,2,3 


dk  dk 

- h  11  :  - 

dt  7  dxj 


-R„ 


dip 

— -  +  ■ 
dx  ,■  dx 


t  v 


K  dk 
a  i.  dx 


~£, 


i  J 


(ID 


and 


ds  ds  s  n  du :  d  f  K  ds 

dt  1  dxj  k  7  dxj  dxj  ^  a e  dxi  y 


-ce2p2/k. 


(12) 


where  K  is  the  turbulent  viscosity  and  cra  ci£,  and  CoWare  numerical  constants. 

For  the  case  study  presented  in  Wang  et  ah,  (6),  the  computational  domain  is  200  m  x  648  m  x 
280  m  with  even  grid  spacing  in  all  directions  (i.e.,  Ax  =Ay  =Az  =  4m).  The  total  number  of  grid 
points  is51xl63x71  and  the  time  step  is  At  =  0.04  s.  Although  not  stated  directly  in  their 
paper,  a  t  =  20  minute  simulation  would  take  approximately  8.3  hours  to  complete.  Figure  4 
shows  an  example  PUMA  model  calculation  of  the  horizontal  wind  field  at  z  =  2m. 
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Figure  4.  A  PUMA  model  calculation  of  the  horizontal  wind  field  at  z  =2m,  i.e.,  wind  velocity 
contours  (a)  and  wind  velocity  vectors  (b)  (from  Wang  et  al.  [<5]). 


2.3  Three-Dimensional  (3-D)  Finite  Control  Volume  Model  (Single  building  and/or  urban 
street  canyon) 

Paterson  and  Apelt  (10)  describe  a  3-D,  flat  terrain,  steady-state,  finite  difference  model  (called 
CITY)  for  a  single  prismatic  building1.  The  CITY  model  time  averages  the  (Navier-Stokes) 
Reynolds  equation  and  the  continuity  equation  to  compute  the  mean  wind  fields.  The  CITY 
model  contains  six  equations  and  the  following  six  unknowns;  the  turbulent  kinetic  energy 
(  k  =  u jU 1 1 2  )  and  its  dissipation  (s),  the  three  mean  velocity  components  (u\  =  u,  ih  =  v,  ih  =  w). 


1  Prismatic  is  defined  as  blocks  with  well-defined  vertical  faces,  where  the  vertical  axes  are  much  longer  than  the  horizontal 
axes  (http://www.onelook.com). 
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and  the  augmented  pressure  (P).  As  presented  by  Paterson  and  Apelt  (10),  the  main  model 
equations  are  as  follows:  Using  regular  summation  notation  for  i,j  =  1,2,3 


dk  _  d 

7  dxj  dXj 


Vt  dk 
crk  dx 


+  v, 


J  J 


dll  j  Su  , 


.  dx  j  dxj 


dll  j 
dxj 


—  £  . 
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and 


ds  _  d 
j  dxj  dxt 


vt  ds 
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dui  i 


.  dx  j  dxj 


diij  s~ 

- ci  — 

dx  j  k 
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and 


duj 
dx ; 
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dxj 


diij 

dx : 

J  y 


dP 

dxj  ’ 


(15) 


duj 

dXj 


0. 


(16) 


Here,  vt  is  the  turbulent  viscosity  and  o/(,  cr.,  c\,  C2,  and  cfl  are  numerical  constants. 

The  above  equations  are  discretized  by  the  finite  control  volume  technique,  i.e.,  partial 
differential  equations  are  integrated  over  appropriate  control  volumes  on  a  staggered  grid  to 
obtain  the  difference  equations.  Here,  hybrid  upwind  differencing  is  used.  The  grid  is  a 
staggered  grid  that  expands  geometrically  away  from  building  faces.  The  method  by  which  the 
pressure  is  calculated  is  known  as  SIMPLE,  i.e.,  Semi-Implicit  Method  for  Pressure  Linked 
Equations  (11).  The  CITY  model  makes  use  of  a  k  -  s  turbulence  scheme  to  resolve  the 
Reynolds  stress.  The  resulting  algebraic  equations  are  solved  by  a  3-D  version  of  the  ADI 
(alternating  direction  implicit)  procedure  in  which  three  sweeps  of  the  solution  domain  (one  in 
each  of  the  coordinate  directions)  are  done  in  each  iteration.  Convergence  takes  about  a  hundred 
iterations  and  requires  about  15  minutes  on  an  IBM  3083E  computer.  Figure  5  shows  an 
example  of  the  wind  velocity  vectors  computed  from  the  CITY  model. 


11 


Figure  5.  Example  wind  velocity  vectors  computed  from  the 
CITY  model  (from  Paterson  and  Apelt  [10]). 


Similarly,  Johnson  et  al.,  (12),  Herbert  et  al.,  (13)  and  Herbert  and  Herbert  (14)  describe  a 
coupled  urban  wind  flow  model  (CITY)  and  a  two-dimensional  (2-D)  thennal  microclimate 
model  (called  SCALAR  and  ENERBAL)  for  city  canyons,  i.e.,  two  buildings.  The  wind  flow 
model  CITY  was  developed  by  Paterson  and  Apelt  (10)  (as  described  above).  For  the  coupled 
urban  model,  the  steady  state  wind  field  is  computed  separately  and  is  maintained  throughout  the 
temperature  simulation.  The  atmospheric  diffusion  equation  in  the  SCALAR  model  as  described 
by  Johnson  et  al.,  (12)  is  as  follows:  Using  regular  summation  notation  for  i  =  1,2,3 


dT  dT 

- b  U ;  - 

dt  dxj 


d 

dx; 


K, 


dT 

dx 


+  S 


i  j 


(17) 


where  T  is  the  air  temperature  at  a  point  in  space,  Kh  is  eddy  diffusivity  for  heat,  and  S  (x,  y,  z,  t) 
is  the  source/sink  tenn.  Here,  only  advection  (by  the  wind  field)  and  diffusion  components  are 
considered. 

In  addition,  the  coupled  urban  model  makes  the  simplifying  assumption  that  the  buildings  on 
each  side  of  the  urban  canyon  are  of  equal  height  and  length,  and  that  the  buildings  are 
rectangular  in  shape,  and  that  each  surface  is  constructed  of  a  homogeneous  material.  The 
SCALAR  and  ENERBAL  models  simulate  the  2-D  temperature  field  within  and  around  an  urban 
canyon  as  a  function  of  the  time  of  day,  time  of  year,  the  wind  field,  location  of  the  city,  the 
canyon  orientation,  the  construction  materials  of  the  buildings  and  street,  and  as  a  result,  the  heat 
fluxes  at  the  building  and  other  surfaces.  Figure  6  shows  a  schematic  of  the  geometry  for  the 
coupled  urban  wind  flow  and  thermal  microclimate  model. 
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Figure  6.  A  schematic  (geometry)  for  the  coupled  wind  flow  and  urban  thermal 
microclimate  model  described  in  Herbert  et  ah,  (13). 


The  coupled  urban  model  domain  is  divided  into  non-overlapping  contiguous  control  volumes 
(e.g.,  240  m  x  632  m  x  100  m).  For  computational  efficiency,  control  volumes  are  selected  to  be 
smaller  close  to  the  ground  and  within  the  canyon  (1-2  m  grid),  and  larger  above  the  buildings 
and  outside  the  canyon  (15-20  m  grid).  The  temperature  in  a  given  volume  of  air  is  treated  as  a 
passive  scalar,  which  does  not  affect  the  wind  flow.  Hence,  the  coupled  urban  model  assumes 
that  the  effect  of  buoyancy  is  negligible  when  compared  to  the  temperature  dispersion  created  by 
the  wind  field.  As  described  by  Herbert  et  ah,  (13),  the  coupled  urban  model  is  implemented  on 
a  Cray  Y-MP8E,  8-processor  super  computer.  On  that  platform,  the  steady-state  wind  field  takes 
about  25  minutes  to  resolve  and  the  combined  thermal  microclimate  model  takes  an  additional 
20  minutes  to  simulate  a  48-h  period.  Figure  7  shows  an  example  result  from  the  coupled  wind 
flow  and  thermal  microclimate  model,  i.e.,  predicted  air  temperatures  (°C)  across  an  urban 
canyon  at  1  pan.  (on  March  15). 
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Figure  7.  Example  results  from  the  coupled  wind  flow  and  thermal  microclimate  model, 
i.e.,  predicted  air  temperatures  (°C)  across  an  urban  canyon  at  1  p.m.  on 
15  March  (from  Herbert  and  Herbert  [ 14 ]). 


2.4  Three  Dimensional  (3-D)  Finite  Volume  -  RANS  (Multiple  building  array) 

Kim  and  Baik  (15)  and  Baik  et  al.,  (16)  describe  their  3-D,  RANS,  CFD  model  with  the 
Renormalization  Group  (RNG)  k  -  £  turbulence  scheme  (17)  to  investigate  non-hydrostatic,  non¬ 
rotating,  incompressible  wind  flows  in  complex  urban  environments.  Their  model  is  based  on 
the  earlier  works  of  Kim  and  Baik  (18)  and  Baik  and  Kim  (19).  The  governing  equation  set  is 
solved  on  a  non-uniform,  staggered  grid  system  (20)  using  a  finite  volume  method  with  the  semi- 
implicit  method  for  pressure-linked  equation  (SIMPLE)  algorithm  (11).  Smaller  grid  sizes  near 
buildings  and  larger  grid  sizes  away  from  buildings  are  used  to  (more  efficiently)  resolve  flow 
and  dispersion  fields.  The  model  equations  as  presented  by  Kim  and  Baik  (15)  are  as  follows: 
Using  regular  summation  notation  for  i,j  =  1,2,3 
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where  the  Reynolds  stress  is  u ( u ,  P*  is  the  deviation  of  pressure  from  its  reference  value,  and 
R  is  an  extra  strain  tenn,  i.e., 


where 


duj  du  j 
ydx  j  dxt 


du,- 


dx  j 


(22) 


(23) 


Here,  ok,  <JE,  C\s,  Cis,  C7„  rjo,  and  f:k\  are  numerical  constants  (18). 

In  the  case  study  described  by  Kim  and  Baik  (15),  a  group  of  buildings  is  embedded  across  101  x 
101  x  41  cells.  The  dimension  of  the  smallest  cell  is  0.3  m  x  0.3  m  x  0.3  m,  which  is  situated  at 
the  edges  of  the  buildings.  The  largest  cell  dimensions  are  1.8  m  x  1.8  m  x  1.8  m.  The  model 
domain  is  63. 1  m  x  63. 1  m  x  28.5  m.  The  time  step  used  in  this  case  is  At  =  0.05s.  The 
computer  model  is  integrated  up  to  t  =  20  minutes  (e.g.,  ~6.7  hours  runtime).  However,  Kim  and 
Baik  (15)  indicate  that  a  quasi-steady  state  in  the  wind  flow  field  is  established  after  t  =  5-7 
minutes.  Figure  8  shows  example  model  results  from  this  calculation. 
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Figure  8.  The  wind  vector  fields  from  the  CFD  model  of  Kim 
and  Baik  (75)  at  (a)  z/H  =  0.5  and  (b)  y/H  =  -0.75  for 
the  case  where  the  initial  wind  direction  is  9=  45°. 


2.5  Three  Dimensional  (3-D)  Finite  Volume  -  LES  (Multiple  buildings;  Central  business 
district) 

Pullen  et  ah,  (21)  and  Boris  (22)  describe  their  3-D,  finite  volume,  CFD  model  (called  FAST3D- 
CT)  with  the  monotone  integrated  large  eddy  subgrid  (MILES)  turbulence  model  (23,24) 
embedded  to  solve  the  high  Reynolds  number,  time-dependent,  Navier-Stokes  equations  for 
mass,  momentum,  potential  temperature,  and  contaminants.  The  time  integration  is  second-order 
accurate  and  has  been  adapted  for  fast  execution  with  very  complex  geometry.  The  CFD 
algorithms  solve  for  slow  but  compressible  flow.  (Note:  While  most  atmospheric  boundary  layer 
models  assume  incompressible  flow,  some  developers  incorporate  compressible  flow  features 
when,  for  example,  detailed  representation  of  thennal  (eddy)  updrafts  and  downdrafts  are  desired 
[25].)  In  addition,  the  model  incorporates  a  complex  finite  volume  algorithm  for  detailed 
building  and  city  aerodynamics.  The  model  physics  implemented  to  compute  the  urban  thermal 
microclimate  includes  solar  heating  of  surfaces  based  on  land-use  data  tables.  The  model 
considers  shadows  from  buildings  and  trees  (depending  on  the  time  of  day)  and  heat  transfer 
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from  building  sides  and  tops  (for  both  daytime  and  nighttime  cases).  Buoyancy  is  included  in 
the  potential  temperature  calculation.  (Note  that  a  list  of  the  governing  equations  for  the 
FAST3D-CT  model  was  not  readily  available). 

Pullen  et  ah,  (21)  presented  contaminant  diffusion  simulations  for  Washington  D.C.  and 
Chicago,  wherein  the  FAST3D-CT  model’s  horizontal  and  vertical  grid  spacing  was  Ax  =  Ay  = 
A z  =  6m.  The  computational  time-step  was  At  =  0.36s.  The  model  grid  embedded  a  1-m 
resolution  building  database.  The  grid  dimensions  were  860  x  580  x  40  levels  for  Washington 
D.C.  and  360  x  360  x  55  levels  for  Chicago.  Typically,  the  FAST3D-CT  model  simulation  of  a 
10  km"  area  at  6  m  resolution  takes  8  hours  on  a  16  processor  super  computer.  Figure  9  shows  a 
contour  plot  of  wind  velocities  for  air  flowing  across  the  Washington,  DC  mall  (from  Boris 
[22]). 


Figure  9.  A  contour  plot  of  wind  velocities  for  air  flowing  across  the  Washington,  DC 
mall  (from  Boris  [22]). 


2.6  Three  Dimensional  (3-D)  Finite  Element  -  RANS  (Single  complex  building 
surrounded  by  a  complex  array  of  trees) 

Calhoun  et  ah,  (26)  present  their  3-D,  RANS,  CFD  model  (called  FEM3)  to  simulate  wind  flow 
and  momentum  around  a  single  complex  building  surrounded  by  a  complex  array  of  trees.  In 
their  study,  numerical  data  are  compared  to  field  measurements.  The  wind  flow  was  assumed  to 
be  neutral,  i.e.,  cloudy,  morning,  or  higher-wind  conditions.  The  turbulence  model  used  is  the 
similarity  /:  turbulence  model,  wherein  the  turbulent  fluxes  are  parameterized  as  proportional  to 
gradients  of  mean  variables.  Canopy  effects  (e.g.,  a  line  of  eucalyptus  trees  to  the  east  of  the 
building)  were  modeled  with  the  addition  of  a  drag  term  in  the  momentum  equations,  following 
Yamada  (27).  The  FEM3  code  uses  a  finite-element  method  (as  discussed  by  Chan  et  al.,  [25]) 
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and  has  been  adapted  for  use  on  massively  parallel  computer  platforms.  The  governing 
equations  for  the  FEM3  model  as  presented  by  Chan  et  ah,  (28)  are  as  follows: 


d 


A  +  plNu  =  -s/p  +  v  . 
dt 


(24) 


V- 


(25) 


_  +  »V6> =  V •  {pCPK°  •  V6>)+  CpN  +  CpA  [k63  ■Vco)-Vd,  (26) 

dt  pC  p  C  p 


—  +  uVa)  =  —  v\pKco-Vco),  (27) 

dt  p  y  ’ 


and 


PM 

RT 


P 


RT 


co 

N 


+  - 


l-a> 


M 


A  J 


(28) 


Here,  u  =  («,  v,  w)  is  the  wind  velocity,  p  is  the  density  of  the  mixture  (e.g.,  dry  air  and  water 

vapor),  p  is  the  pressure  deviation  from  an  isothennal  atmosphere  at  rest  with  corresponding 
density  /?/„  g  is  the  acceleration  due  to  gravity,  0  is  the  potential  temperature  deviation  from 
adiabatic,  co  is  the  mass  fraction  of  the  species  (e.g.,  water  vapor  or  dispersed  contaminant),  K1" 
and  K9,  Ka  are  the  eddy  diffusivity  tensors  for  momentum,  energy,  and  the  dispersed  species, 
respectively,  and  CPN,  CPA,  and  Cp  =  coCPN+  (1  -  co)CPA  are  the  specific  heats  for  the  species,  air, 
and  the  mixture,  respectively.  In  the  equation  of  state  [equation  26],  P  is  the  absolute  pressure,  R 
is  the  universal  gas  constant,  MN  and  MA  are  the  molecular  weights  of  the  species  and  air,  and  T 
is  the  absolute  temperature. 
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Figure  10  shows  the  building  geometry  and  surrounding  array  of  trees  for  the  case  study 
described  by  Calhoun  et  ah,  (26).  Figure  1 1  shows  an  example  result  (modeled  wind  vectors 
versus  experimental  data)  for  the  complex  (building  and  tree)  geometry  shown  in  figure  10. 
Background  momentum  fields  are  also  shown.  Here,  approximately  1.0-2. 5  x  106grid  points 
were  used.  Grid  stretching  allowed  the  finest  grid  spacing  near  the  building  to  be  approximately 
1  m.  The  model  domain  was  400  m  x  400  m  x  80  m.  The  model  simulation  took  approximately 
1  hour  to  complete  on  a  128  processor  super  computer  (i.e.,  the  advanced  simulation  and 
computing  program  [ASCI]  Blue-Pacific  machine). 


Figure  10.  A  schematic  of  the  building  geometry  for  the  case  study  described  by 

Calhoun  et  al.,  (25).  The  circular  and  the  rectangular  shaded  regions  are 
tree  locations  surrounding  the  building. 
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Figure  11.  Example  results  (modeled  wind  vectors  versus  experimental  data)  from  the  FEM3 
code.  White  vectors  are  experimental  data  and  black  vectors  are  model  data. 
Background  shading  represents  modeled  momentum,  where  low  momentum  is 
dark  and  high  momentum  is  lighter  (from  Calhoun  et  ah,  [25]). 


3.  Discussion 


This  section  outlines  some  of  the  main  properties  the  modeling  frameworks  summarized  above, 
to  include  a  discussion  of  difficulties  and/or  deficiencies  associated  with  each  approach. 

3.1  Finite  Difference  Method 

Finite  difference  methods  have  been  around  the  longest  for  numerical  solution  of  partial 
differential  equations  (29).  Some  researchers  consider  finite  difference  methods  to  be  the 
easiest,  most  flexible,  and  most  effective  approach  for  simple  geometries.  Finite  difference 
methods  easily  allow  for  higher-order  schemes  over  regular  grids.  Variable  grid  can  also  be 
implemented  in  a  straightforward  manner  to  allow  for  better  distribution  of  grid  points,  which 
may  help  to  better  resolve  important  atmospheric  scales  and  processes  (30).  In  contrast, 
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numerical  simulation  of  wind  flow  through  complex  geometries,  such  as  those  found  in  urban 
settings,  may  be  quite  difficult  via  finite  difference  methods. 


Nevertheless,  one  begins  by  putting  the  conservation  equation  (e.g.,  mass  conservation, 
advection-diffusion,  etc.)  in  differential  form.  Then,  at  each  grid  point,  the  partial  derivatives  are 
replaced  by  approximations  (e.g.,  forward,  backward,  or  central  differences)  in  terms  of  the 
nodal  value  of  the  functions.  The  result  is  one  algebraic  equation  per  grid  node,  in  which  the 
variable  value  at  that  and  several  neighboring  nodes  appear  as  unknowns  (29). 


As  an  example,  the  conservation  (simplified  Navier-Stokes)  equation  to  describe  the  mean 
concentration  (advection-diffusion)  of  a  scalar  C  can  be  written  as  follows: 

dC  _  dC  d VC 

dt  dx  dz  ’ 


(29) 


where  t  is  the  independent  variable  time,  u  is  the  mean  longitudinal  component  of  the  wind 
velocity,  x  is  range,  z  is  height  above  ground,  and  w'C  is  the  mean  scalar  flux.  The  flux-gradient 
assumption  (31)  suggests  that 

-WC'  =  K  —  ,  (30) 

dz 


Where  K  is  the  scalar  (eddy)  diffusivity.  Combining  equations  29  and  30  yields, 
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In  discretized  form,  this  simple  model  can  be  solved  forward  in  time  using  an  explicit  finite 
differencing  scheme  for  uniform  grid,  i.e., 
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where  i  and  j  are  the  indices  for  the  horizontal  and  vertical  grid,  respectively,  and  n  is  the  time- 
step.  For  an  explicit  scheme,  the  time-step  should  be  small  to  satisfy  the  stability  criterion, 
2KAt/(Az)~  <  1,  as  discussed  by  Press  et  al.,  (32)  (page  838).  Otherwise,  for  larger  time-steps, 
the  numerical  scheme  would  be  unstable  and  the  model  would  not  produce  viable  results. 

3.2  Finite  Volume  Method 

At  the  start,  the  finite  volume  method  uses  the  integral  form  of  the  conservation  equations.  As 
an  example,  the  integral  form  of  a  (simplified)  conservation  equation  to  describe  the  mean 
concentration  of  a  scalar  C  can  be  written  as 

}  CdV  +  \Cu-n  dS  =  \ KVC  ■  n  dS  (32) 

VS  S 


21 


where  J  is  a  volume  integral,  J  is  a  surface  integral,  n  is  the  normal  vector  to  the  surface  of 

v  s 

the  cell,  and  V  is  the  (finite  volume)  divergence  operator  (33).  Here,  finite  volumes  (or  finite 
control  volumes)  are  used  to  discretize  the  equation  set. 

With  the  finite  volume  method  the  model  domain  (space)  is  broken  up  into  a  finite  number  of 
contiguous  volumes  or  cells  and  the  conservation  equations  are  applied  to  each.  Staggered  grid 
is  often  implemented  wherein  cell  centers  are  indexed  j  - 1 ,  j ,  and  j+1  and  cell  edges  are  labeled  j- 
1/2  and  j+1/2.  In  this  manner,  some  variables,  e.g.,  mass  and  energy,  are  evaluated  at  volume 
centers  while  momentum  is  evaluated  at  the  volume  edges.  Interpolation  is  used  to  compute  the 
cell  edge  (surface)  values  in  terms  of  the  cell  center  (nodal)  values.  Surface  and  volume 
integrals  are  approximated  using  suitable  numerical  integration  techniques,  e.g.,  Newton-Cotes 
formulas  (29).  To  its  advantage,  finite  volume  methods  can  be  applied  to  any  type  of  grid,  to 
include  complex  geometries.  The  grid  only  defines  cell  boundaries  and  need  not  be  related  to  a 
structured  coordinate  system.  Nevertheless,  the  selection  of  uniform,  non-uniform,  or  staggered 
grid  will  have  a  significant  effect  on  the  calculation  of  certain  variables.  For  example, 
implementing  a  staggered  grid  may  be  quite  effective  for  pressure  calculations  around  different 
building  geometries,  but  could  create  complications  for  deriving  wind  fields  in  similar 
environments  (25).  Also,  with  a  staggered  grid,  indexing  is  more  complicated. 

3.3  Finite  Element  Method 

The  finite  element  method  is  similar  to  the  finite  volume  method  in  that  the  model  domain  is 
broken  up  into  a  finite  number  of  cells  (which  are  most  often  non-uniform  and  unstructured)  and 
the  integral  conservation  equations  are  applied  to  each.  In  2-D,  the  cells  are  usually  triangles  or 
quadrilaterals,  while  in  3-D  tetrahedral  or  hexahedra  are  often  applied.  The  distinguishing 
feature  of  a  finite  element  model  is  that  the  conservation  equations  are  multiplied  by  a  weighting 
function  before  they  are  integrated  over  the  entire  domain  (29).  An  important  advantage  of  finite 
element  methods  is  the  ability  to  solve  problems  involving  complex  geometries.  Grids  can  be 
easily  refined  by  simply  subdividing  the  individual  elements.  The  disadvantage  of  finite  element 
methods,  as  mentioned  above,  is  that  with  unstructured  grids,  indexing  is  more  complicated. 
Also,  with  finite  element  methods  (computationally)  efficient  numerical  solution  techniques  are 
often  more  difficult  to  find.  Nevertheless,  some  researchers  find  that  finite  element  methods 
provide  a  greater  flexibility  to  model  complex  geometries  than  either  finite  difference  or  finite 
volume  methods  (34). 

Finally,  different  physics/turbulence  algorithms  require  different  kinds  and  amounts  of  computer 
resources.  For  example,  Calhoun  et  ah,  (26)  stated  that  their  RANS  approach  used  about  an 
order  of  magnitude  less  c.p.u  time  than  a  similar  LES  approach.  While  they  found  certain 
advantages  in  using  RANS  over  LES  methods,  nevertheless,  such  models  require  intensive 
computing  capabilities  (even  if  only  for  calculation  of  the  mean  field  variables). 
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4.  Summary  and  Conclusions 


This  paper  presented  a  critical  assessment  of  six  CFD  models.  The  paper  was  derived  from  a 
survey  of  current  numerical  frameworks  to  simulate  wind  flow  and  the  thermal  microclimate 
around  single  and/or  multiple  buildings.  This  study  describes  a  few  complex  CFD  approaches  in 
order  to  better  understand  the  enonnous  task  to  develop  such  software.  CFD  model  summaries 
included  information  on  embedded  physics/turbulence  algorithms,  model  domain,  grid  spacing, 
time-step,  runtime,  and  computer  platform. 

This  paper  showes  that  computer  simulations  of  wind  flow  and  temperatures  in  urban  areas  have 
most  often  been  achieved  via  intricate  and  computationally  intensive  CFD  codes.  The  paper 
provides  an  illustrative  overview  of  current  data  on  this  important  topic.  As  a  result,  it  is 
anticipated  that  the  present  study  will  provide  much  useful  infonnation,  from  which  to  initiate 
several  new  modeling  efforts.  For  example,  it  may  be  advantageous  now  to  explore  alternate 
model  frameworks,  e.g.,  those  that  are  more  computationally  efficient  yet  flexible  enough  for  the 
types  of  future  military  applications  envisioned. 
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